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Microscopic theory of Brownian motion of a particle of mass M in a bath of molecules of mass m ^ 
M is considered beyond lowest order in the mass ratio m/M . The corresponding Langevin equation 
contains nonlinear corrections to the dissipative force, and the generalized Fokker-Planck equation 
involves derivatives of order higher than two. These equations are derived from first principles with 
coefficients expressed in terms of correlation functions of microscopic force on the particle. The 
coefficients are evaluated explicitly for a generalized Rayleigh model with a finite time of molecule- 
particle collisions. In the limit of a low-density bath, we recover the results obtained previously 
for a model with instantaneous binary collisions. In general case, the equations contain additional 
corrections, quadratic in bath density, originating from a finite collision time. These corrections 
survive to order (m/M)^ and are found to make the stationary distribution non-Maxwellian. Some 
relevant numerical simulations are also presented. 

PACS numbers: 05.40.-a, 05.20.-y 

I. INTRODUCTION 

This paper is concerned with the derivation and some properties of the generalized Fokker-Planck equation (FPE) 
for the distribution function f{p,t) of a single stochastic variable p, which differs from the conventional FPE by 
involving p-derivatives of order higher than two, 

n— 1 

Such an equation may appear as a result of a high-order truncation of the Kramers-Moyal expansion of the Markovian 
master equation 

df{p,t) ^ 
dt 



dp'{fip', t)w{p' -^p)- f{p, t)w{p ^ p')} (2) 



n=l 



In general, the naive truncation of this expansion is not a legitimate procedure, but if a problem at hand involves a 
small parameter A, one can approximate the expansion by finite number of terms using an appropriate perturbation 
technique . In case of Brownian motion of a heavy particle of mass M in a thermal bath of light molecules of mass 
m and temperature T, the appropriate small parameter is the mass ratio = mjM. In this case, to order A^ one 
recovers for the particle's momentum p the conventional second-order FPE, 

while going beyond order A^ leads to an equation of the form ([T]) . 

In rare cases when transition rates w in the master equation ([J) are known explicitly, the derivation of the generalized 
FPE HI) is fairly straightforward 0, S, i, i], otherwise it is difficult. A popular approach engages the assumption 
that the fluctuating force in the corresponding Langevin equation is a Gaussian process. In this case perturbation 
analysis is not needed since the terms with derivatives of order higher than two vanish identically and one arrives at 
the conventional FPE which in this case is exact. However, the assumption of Gaussian random force, although 
might seem physically reasonable, should not be taken for granted. In particular, for a Brownian particle it is justified 
only to order A^ , while corrections of higher orders are essentially non-Gaussian and lead to a FPE in the generalized 
form ([1]). This paper is focused on the generalized FPE for the Brownian particle's momentum p to order A^ which 
involves p-derivatives up to order four. 

Despite a few important contributions (see below) , the theory of Brownian motion beyond the lowest approximation 
did not attract much attention in the past, perhaps because of the common belief that higher order corrections are 
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of little importance. However, in recent years the problem has caught some new interest. Beyond lowest order in 
A the Langevin equation for a Brownian particle involves nonlinear dissipative terms, and thus corresponds to the 
description beyond the level of linear response theory. One might hope that the nonlinear Langevin equation and 
the corresponding generalized FPE would enable to capture subtle effects of the interplay of noise and nonlinearity, 
which are completely washed out when one uses the conventional FPE or the corresponding linear Langevin equation. 
This indeed has been demonstrated for a number of problems including Brownian motors 0, 0, BUI and barrier 
crossing [loj . 

The main difficulty related to the generalized FPE and the corresponding nonlinear Langevin equation is that these 
equations usually can not be constructed on a purely phenomenological basis. There is a substantial mathematical 
literature on the Langevin equation with nonlinear dissipation terms. However, for processes with nonlinear dissipation 
thermodynamics provides no hints about the form of fluctuation-dissipation relations, and little progress can be made 
without such relations. Furthermore, an attempt to go beyond the comfortable but artificial assumption of a Gaussian 
random force leaves one, within a phenomenological framework, with no clue how to handle correlations higher than 
of second order. 

These difficulties suggest to derive the generalized FPE from as close to first principles as possible. A successful 
example of such derivation is the van Kampen's method of system size expansion 2J and in particular its application for 
the Rayleigh model of Brownian motion [j, U, @ . In this model a Brownian particle of mass M moves in one dimension 
interacting with the heat bath of ideal gas molecules of mass m <^ M and temperature T through instantaneous binary 
collisions. As was mentioned above, the relevant small parameter is the mass ratio — m/M. Since the momentum 
of the particle P is on average A^^ times larger than that of a bath molecule, it is convenient to work with the 
scaled particles momentum p = XP which is of the same order as the thermal momentum of molecules of the bath 
Pt = \/m kgT. The Rayleigh model is truly Markovian, and the natural starting point is the master equation ^ 
for the distribution function f{p,t), in which transition rates can be readily found explicitly under the assumption 
of binary particle-molecule collisions. Transforming the Kramers-Moyal expansion into the expansion in powers of A, 
one derives to order A^ the conventional FPE 

^ = A^i^./fet), (4) 
where D2 is a second order differential operator 



d 2 9^ } Sv Pt 

op op-' J v27r m 



D2 = la { -^P + Pt^\ ^ To = ^= — , (5) 



u is the number of bath molecules per unit length, pt = y^m/ P is the molecule's thermal momentum, and /3 is the 
inverse temperature f3 — l/ksT. An extension to order A** leads to the generalized FPE in the form 



where the operator D4 involves derivatives up to order four. 



i}fip,t), (6) 



, d 1 a o „ 2 3 ^ 8pl "iPT d'^ , 

D, = j,^--p+ — -p - 2p^— + -—p +——P+——1 (7) 



Note that terms of order A"^ vanish due to symmetry. In what follows I will refer to Eq. ([6]), first obtained in as 
the van Kampen equation. 

One can verify that the Maxwellian distribution fnip) — C exp (— /?p^/2m) is the stationary solution for both the 
standard FPE ([4]) and the van Kampen equation ([6]), 

D2fM{p) - DijMip) = 0. (8) 

However, in contrast to the conventional FPE ([4]), the van Kampen equation ([6]) does not preserve the positivity of 
the solution and therefore can not be an exact equation for any stochastic process. Yet, as an approximation beyond 
lowest order in A, the equation ([6]) is useful and predicts a number of qualitatively new features. For instance, while 
the conventional FPE (U gives for the average momentum {p{t)) — J dpf{p, t)p the closed equation (p) ~ — A^7o (p), 
the van Kampen equation ([6]) predicts the coupling to the third moment (p^): 

j^ip) = -A^o(l - A^) (p) - ^ASoPt' ip')- (9) 
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For initial conditions (p(0)) = and (p^{0)) 7^ 0, the FPE ([4]) gives {p{t)) = for any t > 0, while the equation 
([5]) gives nonzero average momentum for a time interval t ^ Iq^- This prediction was discussed and confirmed by 
numerical simulation inJS]. The van Kampen equation ([6]) has been recently exploited in the context of rectification 
of thermal fluctuations ^ and to study the influence of nonlinear dissipation on the Kramers escape rate fl^ . 

Although proved to be useful, the van Kampen equation ^ is by no means general. It is derived under assumptions 
similar to those for the Boltzmann equation, namely that a characteristic collision time Tc is much shorter than all 
other relevant time scales (which implies small A), and that multiple collisions are negligible (small bath's density). It 
is of interest to derive a generalized FPE from first principles keeping the former assumption, but relaxing the latter. 
Some aspects of this problem were addressed already in pioneering works on the microscopic theory of Brownian 
motion Tl] and further developed in [l^, [3, [13] ■ Perhaps the most elaborate work is the paper by van Kampen 
and Oppcnheim [Tsj . They applied the projection operator technique directly to the Liouville equation for the total 
particle-bath distribution function and derived a generalized FPE of order A**. The coefficients in the equation are 
expressed in terms of rather complicated correlation functions, and no attempt has been made to compare the result 
with the van Kampen equation ([6|) . 

It was noted by Seke that the projection operator method, when applied to the Liouville equation, involves 
some subtlety and may be inconsistent. The alternative way, which we shall follow in this paper, is to apply the 
projector operator technique directly to the equation of motion of the particle, to derive microscopically the Langevin 
equation for the particle's momentum p, and then to apply a standard routine [l| to construct a corresponding FPE 
for the distribution function /(p, t). 

It is generally believed that the two methods of derivation, namely "Liouville — > Fokker- Planck" and "Equation of 
motion Langevin Fokker-Planck" , should give the same result. To lowest order A^ jt is indeed the case: both 
methods lead to the standard second-order FPE ([3]). However, we have found recently .ITj that beyond the lowest 
order predictions of the two methods are different. Moreover, it was found that the second method (from Langevin 
to Fokker-Planck) leads to a generalized FPE with a rather disturbing property, namely, its stationary solution was 
found non-Maxwellian. The intention of this paper is to follow this line in detail to obtain a generalized FPE to order 
A^ in a complete form. The equation we arrive at, namely Eq. (j70p in Section V, involves terms linear and quadratic 
in the density of bath molecules n. For a very diluted bath the latter can be neglected, and the equation is reduced 
to the van Kampen equation ([6]), as expected. The results related to the terms quadratic in n are controversial. 
These terms are absent in the van Kampen equation for a binary collisions model, have the structure different than 
that in the van Kampen-Oppenheim approach, and make the stationary solution inconsistent with Boltzmann-Gibbs 
statistics. Whether these results provide a consistent proof of weak nonergodicity of Brownian motion or are merely 
an indication that a naive perturbation scheme does not apply beyond the Markovian approximation remains an open 
question. As will be discussed, a direct unambiguous verification of non-ergodic effects by numerical simulation may 
be a non-trivial task. 

The possibility of deviations from the Maxwellian distribution has been discussed in literature for a long time. It is 
known that even small deviations may be important for thermally activated processes, in particular for thermonuclear 
reaction rates in astrophysical plasma p^ . Several mechanisms were proposed to justify such deviations, but in our 
opinion none of them are quite satisfactory. For instance, introducing a nonlinear dissipating term into the Langevin 
equation and assuming that the fluctuating force is Gaussian, one generally arrives at a second order FPE with 
a non-Maxwellian stationary solution [l^. However, as we already noted, nonlinear corrections to the dissipation 
force for a Brownian particle is of order higher than A^ . In this case an accurate perturbation procedure leads to a 
generalized FPE of order higher than two, which is a clear indication that the assumption of Gaussian random force 
is not justified beyond the lowest order. 

The plan of the paper is as follows. In sections II and III, the Mazur-Oppenheim version of the projection operator 
technique is applied to derive the nonlinear Langevin equation of order A^. In this part we mostly follow the previous 
paper [20| making some small yet important corrections. The corresponding generalized FPE is constructed in Section 
IV and analyzed in Section V. In Section VI we present the results for the generalized Rayleigh model which allows 
analytical evaluation of all relevant correlation functions. Also, in this section the results of numerical simulation are 
discussed. Summarizing remarks are collected in Section VII. 

II. NON-MARKOVIAN LANGEVIN EQUATION 

Consider a structureless Brownian particle of mass AI immersed in a thermal bath comprised of molecules of mass 
m. It is assumed that the mass ratio A^ = m/M is small and that the bath is initially in equilibrium at temperature 
T. The aim of this and the next sections is to derive the Langevin equation for the particle to order A^. As will 
be shown, such an equation involves a nonlinear correction to the damping term, which is cubic in the particle's 
momentum. 
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The Hamiltonian of the system is 

+ (10) 

= T.ii+ui--x). (11) 

i 

Here x = {xi} and pi are positions and momenta of bath molecules, X and P are those of the Brownian particle, 
Hq is the Hamiltonian of the bath in the field of the Brownian particle fixed at X, and the potential U describes 
the intcrmolccular interaction and the interaction between molecules and the particle. To simplify notations we shall 
consider a one-dimensional problem. The coupling of the particle with hydrodynamic modes of the bath will be 
neglected, in which case the extension to higher dimensions is simple. 

As was already noted, it is convenient to work with the scaled momentum of the particle p — XP, since this quantity 
on average is expected to be of the same order as the typical momentum of a bath molecule pt ~ y/rnksT = ^ m/ (3. 
Writing the Liouville operator C in terms of p has the advantage of extracting the small parameter A explicitly, 

C = £o + A£i, (12) 
m aX op 

Here Fi — ~dU/dxi and F = —dU/dX arc the forces on the i-th bath molecule and on the Brownian particle, 
respectively. 

The operator Cq corresponds to the Hamiltonian Hq and governs dynamics of the bath in the field of the Brownian 
particle which is fixed at X. Since M ^ m one might expect that the force F{t) = e^*F(0) exerted by the bath on 
the particle is close to the force on a fixed particle (pressure), 

F{t) « Fo{t) = e^»*F(0). (15) 

This intuition is implemented in the Mazur-Oppenheim approach plj using the projection operator V which averages 
a dynamical variable A over the canonical distribution p = Z^^exp{—/3Ho), for bath variables sX t — 0, 

VA ={A)= j PA\{ dx, dp,. (16) 

i 

The idea is to decompose the total force F{t) on the particle into a zero centered fluctuating ("random") and a regular 
( "dissipative" ) parts. Using the operator identity 



Jo 

with A = C and B = —VC, one may decompose the force F{t) — e^'^F{0) as follows 

F{t) = F^{t)+ f dT e^^'-'^VCF^ir), (18) 

where F'^{t) = e^^*F and Q = 1 - "P. The factor VCF^t) in the integral in Eq. HH]) can be simplified taking into 
account the orthogonality of V and £o CP^o = 0), and the equality 

(^^FHt)'j=^P{F{0)FHt)), (19) 

which can be verified by integration by parts. As a result, one can write the equation of motion for the scaled 
momentum p{t) = \F{t) in the form [21] 

'^P^^^ = A2 f dr e^(*-^) (-^--p) (^(O)Ft(r)) + XF^t), (20) 



dt 
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where F^t) is a fluctuating force obeying {F^t)) = Ve^^*F{0) = 0. 

The equation (|20p is exact but hardly instructive. To get some progress one needs to expand the fluctuating force 
F^t) = e(^n+^2Ci)t^ in powers of A, 

F^t) = Fo{t) + XF,{t) + F2{t) + 0{X^). (21) 

Using the identity (fTT)) one obtains 

Fo{t) - e^^'FiO), (22) 



/ dhe^o^'-'^^QCiFoih), 
Jo 



F2{t) = f dtie^°^'-'''>QCi f 'dt2e^°'^''~*^-^QCiFoit2). 
Jo Jo 



In contrast to the pressure FqII), the higher order corrections Fi{t) and F2{t) depend on the particle's momentum p. 
Seeking an equation for p(t), one needs to extract this dependence explicitly: 



m 

Here the functions G, are defined as follows 



Fi(t) = ^ / dt'{G^{tX)-{Gi{t,t'))}, (23) 
^^2(i) = J^dt'jyt"[G2{t,t',t")-{G2{t,t',t"))} 
+ 1 fdt' f dt"Go{t^t')\G,{t,t") - {G,{t,t"))}. 



and the operator S is 



Go(i) - ^^o(t), (24) 
Gi(i,ii) = Sit-h)Foih), 

G2{t,ti,t2) = Sit-ti)S{ti-t2)Fo{t2), 



5(t,-tfe) = e^«(*--*'=)^. (25) 



The functions Gi do not depend on the particle momentum p and, as we shall see, all coefficients in the A^-order 
generalized FPE can be expressed in terms of correlation functions of Gq , Gi , and G2 . 

Let us now return to the exact equation of motion Eq. ipO]) and expand the correlation (F(O)F'l'(i)) to order A^, 

{F{0)FHt)) = {F{0)Fo{t)) + X' {F{0)F2{t)) . (26) 

Note that {F{0)Fi{t)) — due to symmetry. Using ([23]) . one can write the above equation in the following form 



{F{0)F^it)) = Go(t) + A2 ( (^)' Gi(0 + - C2{t) 



(27) 



where 



Co{t) = ((Go(O)Go(i))), (28) 

Gi{t)= [ dt'f dt"{{Go{Q)G2{t,t',t"))), 
Jo Jo 

C2(t) - ( dt'l dt"{{Go{0)Go{t-t')Gi{t,t"))). 
Jo Jo 

In the above equations, cumulants {{A1A2 ■ ■ ■ Ak)) are defined in a usual way, i.e. as a part of the correlation 
{A1A2 ■ ■ ■ Ak) which can not be reduced to the product of correlations of lower order. We shall need correlations and 
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cumulants up to order four, which are related to each other as follows: 

(^> - P», (29) 
(^1^2) = (Ai>(A2> + ((AiA2)), 

= (Ai)(A2)(A3) + (Ai)P2A3)) + (A2)P3Ai)) + (A3)P2Ai)) + ((A1A2A3)), 
(^1^2^3A4) = (Ai)(A2)(A3)(A4) 

+ (Al)P2A3A4)) + (A2)PiA3A4)) + (A3)PiA2A4)) + (A4) ((Ai A2 A3)) 
+ ((AiA2))P3A4)) + ((AiA3))((A2A4)) + {{A^Ai)){{A2A^)) + ((Ai A2 A3A4)) . 

The important property of cumulants of functions Gi is that they are linear in the concentration of bath molecules n. 

{{G^GyGu)) -^n. (30) 

This may be demonstrated noticing that Gi is a linear functional of cumulants for the density of bath particles 
N{z,t) — '^iS{z — Zi{t)), where Zi denotes the coordinate-momentum pair (xi,pi) of a bath particle. In turn, one 
can observe that cumulants {{N{zi,ti)N{z2,t2).--N{zk,tk))) of any order k depend linearly on the concentration of 
bath molecules n. For instance, in the expression for the product A^(zi, ti)A^(z27 ^2) one can write the double sum as 
E» J- = J2^^J + T.^=J which gives 

(iV(zi, fi)iV(z2, h)) - {N{ziM)){N{z2,h)) + ^(^(^1 - z,{ti))5{z2 ~ z,ih))). (31) 

i 

Here the second term on the right side is by definition the cumulant {{N(zi, ti)N{z2, ^2))) and is obviously linear in n. 

Substituting the expansion (j27p into the equation of motion (|20p one arrives at a nonlinear and non-Markovian 
Langevin equation of order A^, 

dp{t) 
dt 

Here the memory kernels are 



\' / dr Mi{T)p{t - t) - / dTM2{T)p^{t-T) + \F{t). (32) 



Afi(t) = ^Co(i)-A2^Ci(t) + A2Ac2(i), (33) 
M2{t) = 

and the fluctuating force XF(t) involves three components. 

Fit) = Fo{t) + XFiit) + X''F2{t), (34) 

where Fo{t) is the fluctuating pressure, and Fi{t) and F2{t) are defined by Eqs. (P5|) . The zero-mean term X^F^it) is 
discarded because to order A** it does not contribute to correlations of the fluctuating force XF{t). 



III. MARKOVIAN LANGEVIN EQUATION 

The relative importance of memory effects described by the non-Markovian Langevin equation ([5^ depends on 
how fast the particle's momentum p{t) evolves on the time scale Tc for the decay of the memory kernels Mi(t) and 
M2{t). In what follows we shall assume that the characteristic time Tc does exist. This assumption is not satisfied, 
for example, for the Rubin's model (a heavy impurity embedded in the harmonic chain) where M2(t) is identically 
zero, and Mi{t) decays with time as a power law. It is known, however, that in many cases long-tail effects are indeed 
negligible for sufficiently small A, although the justification may require rather subtle argument [3, [l^l 

The memory kernels Mi(t) are expressed in terms of correlation functions Ci{t) which do not depend on A. Then 
the characteristic decay time of the kernels does not depend on A either, Tc ~ A*'. On the other hand, as follows from 
Eq. (|32p . the characteristic time for relaxation of the particle's momentum Tp is of order A~^ and thus expected to be 
much longer than Tc- This suggests that the non-Markovian equation ([5^ can be expanded in powers of A about its 
Markovian limit. 

One way to make such an expansion is to use in Eq. ([5^ the following substitution [2l[ 



p-{t-r)^p^{t)- J\t'^p^{t'). 



(35) 
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The main contributions to the integrals in Eg. ([32)1 come from the region t < Tc ^ For such r the integral term in 
the right-hand side of Eg. ([55]) is of order p ~ A. Then the nonlinear dissipative term in the Langevin equation ([5^ 
can be written in the local form 

- /" dT M2{T)p^{t - r) = -\^p^{t) f dT Ahir) + 0{X^). (36) 
Jo Jo 

To order A^ this justifies the Markovian ansatz for the nonlinear dissipative term 

- AM dT M2{T)p\t - r) ^ -A'*/(i) / dT M2{t). (37) 



Here the upper integration limit in the right-hand side is taken to infinity since we restrict ourselves to the coarse-grain 
description on the time scale much longer than the characteristic time for bath fluctuations, Tc- 
The same argument for the linear dissipative term gives 



-AM dT Mi{T)p{t-T)^-yp{t) dT Mo{t)+0{\^). (38) 
Jo Jo 

In contrast to the nonlinear dissipation term, the Markovian approximation 

-AM dT Mi{T)p{t - r) -> -X^p{t) / dT Mq{t) (39) 



can be applied for the linear dissipative force only in lowest order A^, in which case one recovers the conventional 
linear Langevin equation 

^^^X^joP{t) + XFo{t), (40) 
with the pressure Fo(t) as a fluctuating force and the dissipation constant 

lo = - r dtCo{t). (41) 
m Jo 

Let us now derive a local form for the linear dissipative term to order A**. Using again the substitution ([35]), one 
can write the linear term as a local expression plus a correction term A{t) 



- Xy dTMi{T)p{t -t) = -X^p{t) / dT Mi{t) + A{t). (42) 

The correction has a form 

A{t)^X'^ I dTMl{T) I dT'p{T'). (43) 



Here the essential integration range is of order Tc A", so that the correction term A ^ X?p ^ X^ . Using the linear 
Langevin equation (|40| . one can write A as follows 

A(i) =-AVp(i) + A^F*(0. (44) 

Here the first term — A'*7*p with 

7*=(£) / rfiCo(t)^ dtCo{t)t (45) 
is a correction to the local linear dissipative force, while the second term 

F*(t)^^ rdTCoiT) f dt'Foit') (46) 
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can be considered as an additional contribution to the fluctuating force. The above expressions for 7* and F*(t) are 
valid for time scale t ^ t^. One can show that in this case F*{t) is a stationary process. Note that the contribution 
F*(t) was overlooked in [20|. 

The result (j44p may be also obtained in a more direct way expanding p{t — r) about t: 

[ dr Mi{t) p{t ^ t) ^ \^p{t) [ dr Mi{t) (47) 



-\^p{t) I dT Mi{t)t + l:\^p{t) [ dr Miiry + ■■■ 
Jo Jo 

To evaluate this expression to order A"* one needs derivatives p*-"-* (t) to order . The first derivative is given by the 
linear Langevin equation (I40p . while for derivatives of higher order Ea. (|40p gives: p^""^' {t)^XF^''-^\t) + OiX^). Then 
the above expansion can be transformed into the form (j42p with the fluctuating term 



X^F^it) = A3 / drA/i(r) \ Fo{t)T - ^Fo{ty + ■■■}. (48) 



Recalling that Mi{t) = ^Co(t) + ©(A^) and noticing that 

1 . 



Fo{t)T-^Foity + ■■■= dt'Fo{t~t')= dt'Foit') (49) 

^' Jo Jt-T 

one recovers F*{t) in the form (|46|) . 

The above results allow us to write the Langevin equation to order A'* in a local form as follows 

^ = -x'j^pit) - xS2PHt) - xm- (50) 

Here the fluctuating force is 

m = Fo{t) + XFi{t) + X^F2{t) + X^F*{t), (51) 

and the dissipation coefficients are 

71 = 70 + A2(57 + A'7^ (52) 



f°° 
70 = - / dtCoit), 
m Jo 



2 3 
h ^ -—J dtC,{t) + £^ / dtC2{t), 
"'JO ^ Jo 

7* = j^dtCoit) j^dtCo{t)t, 

72 = ^ / dtCiit). 



Recall that the correlation functions Ci{t) are defined by Eqs. ([28)) . 



IV. FROM LANGEVIN TO FOKKER-PLANCK 



The aim of this section is to construct a Fokker-Planck equation for the distribution function f{p, t) corresponding to 
the Langevin equation (|50|) . We shall mostly follow the standard procedure [1], but without conventional assumption 
that the fluctuating force is delta-correlated. 

The first step is to assume that /(p, t) obeys a Markovian master equation 



df{p,t) 
dt 



dp'{f{p', t)w{p' ^ p) - f{p, t)w{p ^ p')}- (53) 
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Expressing the transition rates wlpi — > P2) as a function of the initial state pi and the transition length Ap = p2—pi, 
w{pi — *■ P2) — w{pi\Ap), the master equation can be written in the form 

- J d{Ap)[f{p- Ap,t)w{p^ Ap\Ap) - fip,t)w{p\Ap)} 

= j d(Ap)[^{p~Ap,Ap)-^{p,Ap)Y (54) 
where ^(p, Ap) = f{j),t)w{p\Ap). Next, the expansion 



- Ap, Ap) = ^'(p, Ap) + ^ - f -Ap^ j *(p, Ap) 

n— 1 ' ^ 



transforms the master equation into the Kramers-Moyal form 

n— 1 ^ ^ 

with coefficients o;^ given by 



(55) 



(56) 



an{p) = J d{Ap){Aprw{p\Ap). (57) 

Since transition rates w are usually unknown, it is more convenient to work with another representation for q;„. 
Namely, expressing w in terms of the transition probability T{p, t\p\t + t) 

w(p — > p') = lim —T{p,t\p' ,t + t) (58) 

r— ►O T 

one can write the expression ()57p for a„ in the form 

a„(l/) = hm - / dp'ip' - prT{p, t\p',t + r), 

r^O T J 

or 

a„(p) = hm - ([p{t + T)-pit)r). (59) 

r— »0 T \ / 

This expression may be evaluated integrating the Langevin equation for p{t), but this step needs care. Recall that 
the Langevin equation (j50p corresponds to a coarse-grained description with a time resolution much longer than 
correlation time of the random force Tc but much shorter than the characteristic time for the relaxation of the 
particle's momentum Tp. Then in the above expressions the limit t ^ should be understood as Tc <C t ^ Tp. The 

moments (^[p{t + r) — p(t)]"^ in Eq. ([59|) must be first evaluated in the limit r ^ Tc, and only after that the formal 

operation limr^o -^i: ' ') niust be applied: 



an (p) = lim - ( hm ( [p(i + r) - p(t)]" ) \ . (60) 



In what follows, the coarse-grained limit r 3> Tc for brevity will not be indicated. 
Integrating the Langevin equation ([50|) for Tc t Tp 

Pit + t)- Pit) « -[X^jipit) + XS2pHt)]T + A dt'Cit') (61) 



and recalling that ^(t) is a stationary process for t ^ Tc, one obtains from ((60)) : 

ai = -A271P-AS2P', (62) 

a2 = A^limi / dti I dt2(?(ii)?(i2)), 
T Jo Jo 

as = A3 lim - [ dti [ dt2 [ dts (e(ti)e(t2)e(t3)), 
T Jo Jo Jo 

a4 = AMimi / dti [ dt2[ dtsf di4(e(^i)C(^2)^(t3)^(t4)). 

T^O T In Jn . n . n 
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In these expressions the integrals must be taken in the coarse-grained hmit t ^ Tc- 

The next step is to substitute in these expressions the fluctuating force = -fo + + A^F2 + F* retaining 
enough terms to get a„ to order A''. According to ([62|) . the expression for 012 requires the correlation (^^) to order 

A^ 

{ahnt2)) = {Fo{h)Foih)) + X^F,it,)F,{t2)) + X^{Fo{h)F^{t,)) + \^{Foit2)F*{t,)) 

+A2(Fo(ii)F2(t2)> + X''{Fo{t2)F2{h)), (63) 

as requires the correlation (^^^) to order A, 

{a-ti)at2)a-t3)) = X{Foih)Fo{t2)F,{t3))+X{Fo{h)Fi{t2)Fo{h)) + X{F,ih)Fo{t2)Foih)), (64) 

and a4 requires the correlation (^^^ to order A°, 

{ahmt2mhmh)) = {Foih)Fo{t2)Foih)Fo{U)). (65) 



To extract the dependence on p one has to express Fi and F2 in terms of p-independent functions Gi (t) , see Eq. (|23p . 
Then the correlation functions take the forms: 

mimt2)) ^Co{tl,t2)+X^m-^{C2{h,t2) + C2{t2,h)} (66) 

+X^ {Ci{ti,t2) + Ciit2,h) + C3it2,h)}+X''^{Ce{ti,t2)+Ceit2,ti)Y 

(e(il)e(i2)e(i3)) - A-|C4(tl,t2,t3) +C4(tl,t3,t2) +C4(i3,t2,tl)|, 

{^{hMhMtsMh)) ^ Co{h,t2)Co{t3,h) + Co{tut3)Co{t2,h) + Co{h,U)Co{t2,t3) + C5{ti,t2,t3,h). 

Here the functions Ci are defined as follows 

Co{ti,t2)^{{Go{ti)Go{t2))), (67) 

Ci(ti,<2)- f'dt' f dt"{(Go{ti)G2{t2,t',t"))), 
Jo Ja 

C2{tl,t2)^ r'dt' [ dt"{{Go{tl)Go{t2~t')Gi{t2,t"))), 
JO JQ 

Cz{tiM) = / 'dt' f 'dt"{{Gi{ti,t')Gi{t2,t"))), 
Jo Jo 

C^itiMM) = / 'dt'((Go(ii)Go(t2)Gi(t3,t')», 
Jo 

C^itiMMM) = ((Go(ii)Go(t2)Go(t3)Go(t4)», 
C&{hM)^ f'dt'Coit') r dt"Co{ti,t"). 

Jo Jt2-t' 

Note that the two variable function Go(ii,t2) = ((Go(ii)Go(i2))) defined above, and the single variable function 
Co{t) — ((Go(O)Go(i))) which we used in the previous sections, are related as Go(ti,t2) = Go(|ti — t2\)- 
Substitution of correlations (|66p into leads to the following results 

'72 p', (68) 

' 74 + A —75 + A —76, 
m m 



ai 


= -A^7ip- 


Oi2 


= aS + a' 




= A —77, 




TO 


a4 
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Here the dissipation coefficients 71 and 72, are defined above by equations (l52l) . 



73=limi / dhl dt2Co{ti,t2), (69) 



r^O r 



JO 



74=limi f dtJ dt2\2C\{ti,t2) + C3{tut2)}, 
T Jo Jo ^ J 

75 = lim- / dtJ dt2 2C2{ti,t2), 

r Jo Jo 

76 = lim i / dtif dt22Co{h,t2), 

T Jo Jo 

77 = lim - / dtij dt2 dt33C4{ti,t2,t3), 

T Jo Jo Jo 

1 

78=Hni- / dti dt2l dtsj dtiC^{ti,t2,h,U). 



"'0 



Recall again that in these formulas the integrals must be evaluated in the coarse-grained limit t ^ Tc. 

Note that in the formula for 73 we have discarded the terms involving the products of cumulants, like 
{{Fo{ti)Fo{t2))){{Fo{t3)Fo{t4))). Such products depend on two time differences, and the corresponding contributions 
to the four-dimensional time integral in the expression for are quadratic in r. They therefore vanish when the 
operation linv^o 7(- ■ ■ ) is applied. 

V. GENERALIZED FOKKER-PLANCK EQUATION 



Substituting the results (|68p for a„ into the Kramers-Moyal expansion ([SG]) one arrives at the Fokker-Planck equation 
of order in the following form 

^ {X'D2 + X'D, + X^D*2}fiP,t). (70) 

Here the differential operators D2 and D4 have the same structure as for the van Kampen equation 

d 52 

D2=a,-p + a2^, (71) 

d d 92 (9^ 

Di = bi — p + b2j-p^ + b:i— + bi—p' + h—p + hoTr^. (72) 
op op Op^ Op'' Op'^ op^ 

The difference with the van Kampen equation is the presence of the operator L'Jj originating from non-Markovian 
corrections 7* and F* . It has the same structure as D2 

+ (73) 

but as we shall see, scales differently with the bath density n. 
For the operator D2 the coefficients are ai = 70 and 02 — 5 73, 

ai = ^ / dtCoit), (74) 
■m Jo 

02 = TT i dtij dt2Co{ti,t2). 



2t Jo 







The expression for 02 may be simplified recalling that Co{ti,t2) — Co{\ti — t2|) and using the coarse-grained limit 



02- - / di2 / dhCo{tl,t2) = - f dt2{T-t2)Co{t2)^ I 

T Jo Jo ''' Jo Jo 



rt2 

dtCo{t). (75) 
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This gives the conventional relation 

^ = ^, (76) 
a2 rn 

which guarantees that the Maxwellian distribution fM{p) ~ C exp{—l3p^ /2m) is stationary for the operator D2: 

D2fM{p)^0- 

For the operator D4 the cocfBcients are 61 — Sj, 62 — 72, 63 = Ji/^rn^, = 75/2771, 65 = 77/31771, and bg = 78/4! 
Using the results ([52|) and (l69|) for 7^, one can write the coefficients bi in terms of correlation functions as follows: 

dtCiit) + ^ f dtC2{t), (77) 





2 




771^ 


&2 = 


JL 








1 


&3 = 




277l' 




1 


64 = 






771 




1 


65 = 


2m 






&6 = 


4! T 



lini- / dti/ dt2l dt3C4{ti,t2,h), 
" ~ 3 "'0 "'0 

dti dt2 dts dt4C5{ti,t2,t3,t4). 
Jo Jo Jo 

Here the functions Ci are defined by relations ([28|) and (l67|) , and the integrals must be evaluated in the coarse-grained 
limit T ^ Tc. 

We did not attempt to give a general prove that D4, with coefficients given above, satisfies the stationary relation 
-C'4/m(p) = 0. Instead, in the next section we evaluate coefficients bi explicitly for the exactly solvable generalized 
Rayleigh model. In this case D4 is found to be the same as for the van Kampen equation ([6]), and therefore the 
relation D4fi\i{p) = is satisfied. 

Consider at last the coefficients for the operator • 



It is proved in the Appendix that 



ci=7*=(^^j dtCoit) dtCo{t)t, (78) 
C2 = — 76 = — lim - / dti dt2C6iti,t2). 

2m m T^O T Jn Jn 



C2 2m 



Because of the factor two in the denominator the Maxwellian distribution /jv/ ~ exp{—f3p'^ /2m) is not stationary 
for the operator I?2, D2fM{p) 7^ 0. The validity of the relation ([75)) can be also verified directly, for instance, for 
Co{t) ^ exp{—t/Tc). One might hope that keeping terms of even higher orders would restore Maxwellian equilibrium. 
However, the distribution /a/ (p) does not depend on A, which suggests that it must satisfy each term of the A-expansion 
separately. 

Note that coefficients ci and C2 are quadratic in cumulants Co{t) = {{Fo(ti) Fo{t2))) , and therefore quadratic in the 
bath density, D2 ^ . In contrast, the operators D2 and D4 are linear in n. Therefore, in the low-density limit D2 
may be neglected, and the generalized FPE ([70]) is reduced to the van Kampen equation ([6]). On the other hand, one 
may expect that for sufficiently high density of the bath the operator D4 ^ n may be neglected compared to ~ n^, 
which results in a conventional FPE 

9f(P^'^-X^lA,lp + A2^\fip,t), (80) 



dt \ dp^ ' " dp^ J 

but with modified coefficients 

Ai=ai + \^ci,A2^a2+\^C2. (81) 
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The stationary solution of this equation is the Maxwelhan distribution exp [— /3*p^/2m] with the inverse temperature 

/3 -m— — (82) 

which is smaller than that for the bath, /3. Indeed, recalling that 01/02 = P/m and Ci/c2 = [3 /2m, one obtains to 
order 

^ = l-X^^^l-\^^ r dtCo{t)t. (83) 
P ai m Jq 

Needless to say, the prediction that the temperature of a Brownian particle is higher than the temperature of the 
bath is in contradiction with basic assumptions of equilibrium statistical physics and must be subjected to a thorough 
scrutiny. It might be instructive to consider a specific model. 



VI. GENERALIZED RAYLEIGH MODEL 

For the generalized Rayleigh model [l^l it is possible to evaluate the coefhcients in the nonlinear Langevin equation 
((50)1 and generalized Fokker-Planck equation ([70)l analytically. In this model the bath molecules do not interact 
with each other, while the Brownian particle interact with molecules not through instantaneous collisions, as in the 
original Rayleigh model [3) S but via a continuous parabolic repulsive potential. Namely, when the distance 
between a molecule and the particle |a;i — X| is larger than a given length R the molecule moves freely. But when the 
molecule enters the "interaction zone" \xi — X| < i?, it experiences a repulsive parabolic potential ■^k{xi — X)'^, where 
X = X — R for a molecule approaching the particle from the left, and X = X + R for a molecule from the right. 

As was shown in for this model all relevant correlation functions can be evaluated exactly. In particular, the 
correlation functions Ci{t), Eqs. (j28p . which determine dissipative coefficients in the nonlinear Langevin equation 
(ISO]) , read 

Co{t) = LJvp^m^^CoiuJt), (84) 
Ci{t) = uj v pt m ^i{ujt), 
C2{t) = -ujvp% £,2{ujt). 



Here u = yjk/m is the inverse collision time, v is the number of molecules per unit length, pT = \/ m/ (3 is the thermal 
momentum of a molecule, and dimcnsionlcss functions ^i{x) are 

2 

^o(a;) — — 6{tt ~ x){sin x + {-n — x) cos a:}, (85) 



S,i{x) — 9{t: — x) siv? x, 



/2tt 

£.2{x) — 9{tt — .T){sin^ X + x{tt — x) sinx}. 



/2tt 

where 9{x) is the Heaviside unit step function. Substitution of these results into Eqs. gives the dissipative 

constants for the nonlinear Langevin equation (j50p 

^ = X'j,p{t) - XS2pHt) - Xm (86) 

in the form 

71 = 7o + A2^7 + AV, (87) 



70 = / dxiQ{x) 



m Jq V27r m 

= -'^ dx[2^,{x)+U^)] 
m Jo 



8 vpx 



/2n m 



^2^2 roo ^00 ^2^2 

7* = 1^ / dxS,fi{x) / dx^o{x)x = S 1, 

^m'^ Jo Jq ^rn^ 

72 = — ^ / dx^i{x) 



mpT Jo 3%/27r mpx' 
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Consider now the generalized FPE ([70)1 

df{p,t) 
dt 



[x'D2 + X^Di + X''D*yf{p,t). 



The coefficients in the operators Di are determined by correlations (j67p . They can be evaluated in the same manner as 
the correlations ([M|) . Then one can show that the operators D2 and D4 coincide with those for the original Rayleigh 
model with instantaneous binary collisions and are given by Eqs. ([5]) and ([7]), respectively. 

The operator Z?^ originates from non-Markovian corrections and does not appear in the generalized FPE for the 
original Rayleigh model with instantaneous collisions. It has the form 



with coefficients 



ci=7^ = 8^, C2 = 16^. (90) 



In accord with the general prediction (|79p . one observes that Ci/c2 = /3/2m, and therefore the Maxwellian distribution 
/m(p) — Cexp(— /3p^/2m) is not stationary for the operator, I?5/m(p) 7^ 0. 

Recall that D4 and D2 scale differently with the molecular density i': D4 ~ i', D2 ^ . This suggests that for 
sufficiently high density, the operator D4 can be dropped. Inspecting Eqs.© and ((90)) . one finds that the ratio of 
terms generated by to those produced by D4 is of order 

iV=^. (91) 

This is an important parameter of the problem and has a meaning of average number of molecules simultaneously 
interacting with the particle. For the case <C 1, corresponding to the limit of binary collision, the operator D\ 
can be neglected, and the generalized FPE is reduced to the van Kampen equation The opposite limit S> 1 
corresponds to multiple collisions. In this case the operator D\ is expected to be more important than D4, and the 
generalized FPE to be reduced to the form ([80|l . As discussed in the previous section, this equation has a Maxwellian 
solution with the inverse temperature /3* smaller than that for the bath /3. The ratio /3*//3 is given by Eg. ([55)1 . which 
for the given model takes the form 

^ = 1 - ^N)?. (92) 

We attempted to verify this prediction in a numerical experiment using a simulation setup similar to that described 
in detail in [9]. In our simulation two Poissonian sources of in-going molecules are located sufficiently far from the 
particle to mimic an infinite thermal bath with a Maxwellian velocity distribution of molecules. The particle's position 
is initially fixed and the bath is allowed to equilibrate in the field of the fixed particle for a time t^q- Then, at the 
moment t = the particle is released with an initial velocity 1^(0), and the relaxation of the mean-square velocity 
(F'^(t)) towards its stationary value is calculated. Fig. 1 shows the result for the mass ratio = 10""^ and the 
density of the bath corresponding to TV = 100. The velocity is in units of the thermal velocity of bath molecules 
vt = l/\/rnP, so that according to equilibrium statistical mechanics (y^(<)) should approach the equilibrium value 
A-^ = 0.001. Instead, the stationary mean-square velocity of the particle appears to be shghtly higher. However, the 
exceeding is less than 0.1 percent which is orders of magnitude lower than predicted by Eg. ([5^ . The smallness of the 
alleged deviation from the Maxwellian average requires very large number of sampling runs, which in our case was 
about 3 • 10^. Also, the result seems to be sensitive to the value of the integration time step Ai. We used the velocity 
Verlet algorithm with At ^ 0.001 in units of the collision time cj~^. Simulations with larger At show no signs of 
deviations from Maxwellian statistics. 



VII. CONCLUDING REMARKS 



In this paper the Fokker-Planck equation for the Brownian particle is attempted to be derived beyond the lowest 
order in the mass ratio A^ = m/M . Unlike the approach of van Kampen and Oppenheim p^ , the projection operator 
technique was applied to the equation of motion for the particle, rather than to the Liouville equation. The results 
seem both encouraging and controversial. 
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2 4 6 8 10 12 14 

Time 

FIG. 1: The relaxation of the mean-square velocity of a Brownian particle {V'^{t)) for the generalized Rayleigh model. 
Parameters of simulations are = 0.001, N = 100, teg = 7.5. Time is in units of the collision time uj ^, and velocity is in 
units of the thermal velocity of bath molecules vt- Zoomed data shown in the inset suggest that the stationary value exceeds 
the Majcwellian average (1/^) = = 0.001. 

On the one hand, the van Kampen equation, originally designed for the special case of instantaneous collisions, is 
recovered from first principles for a low-density bath. Therefore, in the Markovian limit the method is correct. 

On the other hand, it is found that in general case the equation contains the additional operator A^DJ originating 
from non-Markovian corrections which are inevitable for any model with finite collision time. These corrections are 
found to make the stationary solution non-Maxwellian, D^ fMiv) ^ 0, which is, of course, a very disturbing result. 
Brownian motion of a massive particle coupled to an infinite bath in thermal equilibrium, with the Maxwellian 
velocity distribution for bath molecules, is often considered as a classical example of an ergodic process: all accessible 
microstates are equally probable over a long period of time, and therefore the equilibrium state is the Gibbs canonical 
ensemble. To order this anticipation is supported by solutions of both Langevin and Fokker-Planck equations. 
To order A"', the original Rayleigh model with instantaneous collisions leads to the van Kampen equation which also 
has the Maxwellian stationary solution. Note that in this last case the ergodic behavior can not be deduced from 
the Central Limit Theorem or the assumption of Gaussian noise imposed on the Langevin equation (which predict 
the FPE of second order). Instead, one can show that the van Kampen equation is equivalent to a master equation 
with transition rates obeying the detailed balance condition As well known, in this case the applicability of 
Boltzmann-Gibbs statistics can be readily proved 0. However, whether or not the detailed balance holds in general 
case is not known, and a truly dynamical justification of Boltzmann-Gibbs statistics still seems to be missing [23]. A 
number of generalizations of Boltzmann-Gibbs statistics has been discussed in recent years, including those derived 
directly from underlying dynamics (see for instance [IJiHEIl^ and references therein). However, these generalizations 
usually imply an essential departure from basic assumptions of Boltzmann-Gibbs statistical mechanics such as inelastic 
collisions, nonextensivity. Levy statistics for the bath, etc. In contrast, the result Dl^fMiv) 7^ originates merely from 
finite duration of collisions, which suggests deviations from Boltzmann-Gibbs statistics under a much wider range of 
conditions. 

This radical prediction is in contradiction with the results of van Kampen and Oppenheim 15]. They obtained a 
generalized FPE in which the operator (in our notations) has the form of D2 squared, and therefore the stationary 
solution is Maxwellian. This casts some scepticism about the ability of the perturbation approach exploited in this 
paper to treat non-Markovian effects. The method is based on the seemingly straightforward assumption of wide 
separation of time scales for the particle's momentum Tp ~ A~^ and for correlation functions of the fluctuating force 
Tc ~ A*^. However, for a finite A, the relation Tp ^ is not necessarily true, even if the coupling with slow bath's 
collective modes is negligible [IJ]. On the other hand, the van Kampen-Oppenheim approach also relies on the 
assumption ^ Tc, and it is not clear why the two approaches give different results. Numerical simulations of the 
generalized Rayleigh model apparently suggest some deviation from the Maxwellian statistics, but much smaller than 
the theory predicts. It should be stressed, however, that the presented theory is asymptotic and implies the weak 
coupling limit, which is difficult to approach in a numerical experiment. One may hope that further, more accurate 
numerical modeling would shed some light on these questions. 
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Appendix 

In this Appendix we prove the relation ([79|l : ci/c2 = (3/2m. Using the definition of ci and C2, Eg. ([78|) . the relation 
can be written in the form 



f dti / dt2Ce{tiM) — > 2t dtCoi 
Jo ^=^^<= Jo 



OO 



(t) dttCoit). (Al) 







or taking the derivative 

dtC6{t,T)+ I dtCe{T,t) — > 2 / dtCo{t) I dttCoit). (A2) 



Asymptotic relations (|Aip and (|A2[) are equivalent but the latter is easier to prove. Let us evaluate two integrals in 
the left hand side of the relation (jA2p . 

Consider the first integral Ii — dtC%{t,T). Recalling the definition of C% 



Ceih,h) = r dt' Co{t') f ' dt" Coilh - t"\) 

Jo Jt2-t' 



one can write 



h = dt dt'Co{t') dt" Ca{\t-t"\) 

Jo Jo Jr-t' 



dt'Coit')f dt"U dtCait" -t)+ I dtCa{t-t")} 



Jr-t' ^JO 

dt'Co{t')f dt" f dt"'Co{t"')+ f dt'Coit') f dt" f dt"'Co{t"'). 

Jr-t' Jo Jo Jr-t' Jo 

Introducing the function 

$(0 = / dt'Coit') 
Jo 

the above expression can be written as follows 

h= f dt'Coit') f dt"<^it")+ I dt'Coit')! 

Jo Jr-t' Jo Jo 

Noticing that ^'(i) = C'o(i), it is convenient to integrate the last term by parts: 

I dt'Coit')! dt" ^it") + <^iT) I dt^it)- [ dt^'^it). (A3) 

Jo Jr-t' Jo Jo 



Consider now the second integral in the left hand side of (|A2|) 

/2= / dtCeiT,t)= f dt f dt'Coit') f dt"Coi\T-t"\). 
Jo Jo Jo Jt-t' 

Since t" < t < t the symbol of absolute value may be omitted, and a change of variables gives 

ct" f-r J rr ft" 



h ^ I dt' I dt" Coit" - t') / dt'"Coit"') = - / \t' t'-^ f dt" Coit" - t') I dt'" Coit" 
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Evaluation of the derivatives gives 

/2 = Tdi'i' Co{t - t') fdt" W)+ fdt' t' Co{t') $(r - t')- fdt' t' fdt" Co{t" - t')Coit"). 
Jo Jt' Jo Jo Jt' 

Integrating by parts, the last term in this expression can be presented in the form 



dt't' dt" Co{t" ^t')Co{t") = -^{t) / dt'^{t')+ / dt^\t), 
Jt' JO JO 



so finally for I2 one obtains 



l2= / dt't' Co{T-t') dt"Co{t")+ / dttCo{t)^{T-t)~<^>{T) / dt^t)+ / dt<i>^{t). 
Jo Jt' Jo Jo Jo 

Together with (jA3|) this gives for the left hand side of relation (|A2|) 

/•T /"T /"T /"T /"T 

h+h^ / dt't' Co{T~t') dt"Co{t")+ / dt'Coit') dt"<i>{t")+ / Co(t) $(r - t). 

Jo Jt' Jo Jr-t' Jq 

In the limit t ^ Tc the first term in this expression vanishes, while the second and the third terms both equal 
$(00) f^dttCoit): 

h+h — > 2$(oo) / dttCoit). (A4) 
But this is just the right hand side of the relation (jA2p . which thus is proved. 
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